% 蹦床疲劳建模 - 简化模型
clc; clear;

%% 参数设置
g = 9.81;                          % 重力加速度
Nx = 5; Ny = 5;                    % 网格尺寸
N_nodes = Nx * Ny;                % 总节点数
T_total = 10;                     % 模拟总时间（s）
dt = 0.01;                        % 时间步长
Nt = T_total / dt;                % 时间步数

% 节点参数（均匀）
m_node = 1.0;                     % 单个节点质量（kg）
k_node = 5000;                    % 刚度（N/m）
c_node = 20;                      % 阻尼（N·s/m）

% 初始化变量
z = zeros(N_nodes, Nt);          % 位移
v = zeros(N_nodes, Nt);          % 速度
a = zeros(N_nodes, Nt);          % 加速度
F = zeros(N_nodes, Nt);          % 外力

%% 运动员参数
athletes = [
    1, 68, 2, 2, 0.5;  % [编号，质量kg，x节点，y节点，起跳时刻]
    2, 55, 4, 2, 2.0;
    3, 72, 3, 3, 4.0
];
N_ath = size(athletes,1);

% 起跳参数
F_peak = 4 * g;                  % 峰值为体重的4倍
T_contact = 0.4;                 % 蹬压时间
omega = pi / T_contact;

for i = 1:N_ath
    m_i = athletes(i,2);
    x = athletes(i,3);
    y = athletes(i,4);
    t_start = athletes(i,5);
    
    node_idx = sub2ind([Nx,Ny], x, y);  % 位置转索引
    t_idx = round(t_start/dt):round((t_start+T_contact)/dt);
    
    % 外力模型：正弦型蹬压力
    for j = 1:length(t_idx)
        t_j = (j-1)*dt;
        F(node_idx, t_idx(j)) = F(node_idx, t_idx(j)) + m_i*g + m_i*F_peak*sin(omega*t_j);
    end
end

%% 模拟蹦床动力学
for t = 2:Nt
    for i = 1:N_nodes
        a(i,t) = (F(i,t) - c_node*v(i,t-1) - k_node*z(i,t-1)) / m_node;
        v(i,t) = v(i,t-1) + a(i,t)*dt;
        z(i,t) = z(i,t-1) + v(i,t)*dt;
    end
end

%% 疲劳损伤计算（基于最大应力）
% 应力映射：σ = k*z（简化）
sigma_max = max(abs(k_node * z), [], 2);    % 每个节点最大应力
a_SN = 12; b_SN = 3;                         % 假定 S-N 曲线：logN = a - b*logσ
N_cycles = 10;                              % 假定每节点受10次重复载荷

D = zeros(N_nodes,1);                       % Miner损伤值
for i = 1:N_nodes
    if sigma_max(i) > 0
        N_i = 10^(a_SN - b_SN*log10(sigma_max(i)));
        D(i) = N_cycles / N_i;
    end
end

%% 可视化损伤分布
D_grid = reshape(D, [Nx, Ny]);
figure;
imagesc(D_grid);
colorbar;
title('蹦床疲劳损伤热图（Miner D值）');
xlabel('X'); ylabel('Y');
